Disorder and Interactions in ID Systems 
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We present a new numerical approach to the study of disorder and interactions in quasi- ID systems 
which combines aspects of the transfer matrix method and the density matrix renormalization group 
which have been successfully applied to disorder and interacting problems respectively. The method 
is applied to spinless fermions in ID and the existence of a conducting state is demonstrated in the 
presence of attractive interactions. 
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I. INTRODUCTION 

It is well established that in the presence of disor- 
der electron wavefunctions can become localized. Con- 
siderable numerical work has been carried out for non- 
interacting systems with results reaching a reasonable 
consensus: theory and experiment are in general qual- 
itative agreement. However, in 3D the calculated value 
of the universal critical exponent is markedly larger than 
the empirically measured valuei. This seems to sug- 
gest that an essential factor is missing from calcula- 
tions: the obvious candidate is the electron-electron in- 
teraction. Furthermore, some have claimed to observe a 
Metal-Insulator transition in 2D contrary to the widely 
accepted scaling theory of Anderson localization^. This 
is often accredited to the effect of interactions. Hence 
during the last 10 years attention has been switching to 
this more difficult case. The central problem is that the 
model becomes a many-body system and so the Hilbert 
space grows quickly with system size. This renders an ex- 
act numerical calculation far beyond computational ca- 
pabilities. Nevertheless, several studies have been accom- 
plished, these suggest inclusion of interactions may yield 
non-trivial behavior. 

Shepelyanskja* performed calculations on two interact- 
ing particles. In ID, interactions caused a large en- 
hancement of localization length. Other work showed 
that in 2D the effect is possibly stronger leading to 
delocalizationi. However, some caution is required as the 
method fails to reproduce known non-interacting results 
when interactions are switched off. 

The most successful method for treating the finite 
density problem is the Density Matrix Renormalization 
Group (DMRG) approach^. This works by performing 
a direct diagonalization but reducing the Hilbert space 
by systematically discarding basis states that do not con- 
tribute significantly to the ground state. Applying this 
method to the Anderson interacting model (defined in 
equation a delocalizcd regime was found for attrac- 
tive interactions^. In more recent papers by the same 
authors, it was noted that interesting physics is washed 
out in the averaging process. Charge reorganizations can 
be seen as electrons on a chain shift from the Mott insu- 
lator limit (strong interactions) to the Anderson insula- 



tor limit (strong disorder)^.. Extensions of DMRG to 2D 
have encountered difficulties. 

We have developed a new method incorporating some 
of the ideas of DMRG and the transfer matrix method 
successfully used in the non-interacting cas oV . Sec- 
tion [H] describes the method and section II I II discusses 
the application to a model of spinless fermions. 



II. THE NEW METHOD 

Like DMRG our approach is based on a tight-binding 
method and has a similar potential usefulness. It can be 
readily applied to describe any ID or quasi-lD system, 
provided interactions are nearest-neighbor. 

A. The Hamiltonian 

As this is a many-body problem it is natural to work 
within the second quantization formalism of quantum 
mechanics. This allows the Hamiltonian to be written 
in terms of particle creation cf and annihilation c i oper- 
ators for site i: 

i i 

i i 

The first two terms constitute the standard Anderson 
modelii used widely in the study of disorder induced lo- 
calization. The additional U term represents the nearest- 
neighbor interaction. If neighboring sites are occupied 
then the two particles experience a repulsive force (as for 
electrons) or possibly attractive force of strength U . Set- 
ting U = reverts the system to the many independent- 
body situation (i.e. the non- interacting case). 

The final term represents the chemical potential /i. It is 
necessary as this method works within the grand canon- 
ical scheme in which a range of particle numbers will be 
considered. The value of the parameter fi corresponds to 
the energy required to add a particle to the system, and 
thus controls the particle density of the system ground 
state. As with most numerical studies of Anderson local- 
ization zero temperature will be assumed. 
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B. The Recursive Method 

Our method tackles the problem of the exponentially 
growing Hilbert space by reducing the number of basis 
states, restricting the focus to the ground state. This 
works in conjunction with a recursive procedure that ex- 
tends the chain by successively adding new sites. Open 
boundary conditions must be used. 

oooooooo 

FIG. 1: The recursive procedure adds new sites to both ends 
of a ID chain at each iteration. 

For each iteration: 

• A site is added to each end of the chain (fig. 
and basis states are constructed. At first sight it 
may appear simpler to add a new site to one end 
only. However, it turns out that for the purposes 
of measuring the degree of localization it is much 
more natural to add sites to both ends of the chain 
in the same iteration (section III C() . 

• For each particle number with remaining basis 
states a Hamiltonian matrix is found. After the 
matrix elements have been calculated, the Hamil- 
tonian is solved and the desired quantities are ex- 
tracted. 

• Finally, a proportion of the resulting eigenstates 
are thrown away according to some criterion. The 
remaining states are used to form the basis at the 
next iteration - a chain with two more sites. 

There is no fundamental reason why this process can- 
not be repeated to very large chain lengths. The following 
sections detail the mathematics of this procedure. 

1. Expressing the Basis States 

For a one-dimensional chain, with L sites and one elec- 
tron per site, there are 2 L basis states. In order to reduce 
the Hilbert space, this new method relies on the fact that 
it is possible to express the states for a chain of length 
L in terms of the energy eigenstates |$ L_2 ) of the chain 
of length L — 2 (i.e. the same chain without the two end 
sites). Formally, the L site Hilbert space is a product of 
the Hilbert space for the L — 2 site chain with the vector 
space associated with the two new sites. 

The eigenstates |$ L_2 ) can be used as an orthogonal 
basis for the inner Hilbert space because they are eigen- 
vectors of the previous iteration Hamiltonian such that 

H L - 2 \$f- 2 ) = Ei\$f- 2 ). (2) 



The outer Hilbert space, associated with the two new 
sites, is spanned by four basis states. This is readily 
seen by considering particle number occupancy represen- 
tation: |0 - - - 0>, 1 1 - • - 0> , |0 - - - 1) and Thus for 
every eigenstate |<Ff ~ 2 ) of the L — 2 site chain, there are 
four corresponding basis states for the L chain: |0$f _2 0), 
|l$f" 2 0), |0$f" 2 l) and |l$f~ 2 l). 

Consequently, a general state for the L chain with TV 
electrons, \^^' N ), may be written as a linear combination 
of basis states in the following manner: 

l*£' w >= EMosf-^o) 

i 

3 

+ Y,d nk \l^ L k -' 2 ' N - 2 l). (3) 

k 

In fact, as this equation indicates, it is only necessary to 
consider the subset of |$ i_2 ) states which, when com- 
bined with two new end sites, have a total of N electrons. 
The reason is that particle number is a good quantum 
number for this Hamiltonian. 



2. Calculating the Hamiltonian Matrix 

Thus basis states can be grouped according to particle 
number and a separate Hamiltonian can be calculated for 
each. This can only be accomplished by first expanding 
each of the four types of basis states for N particles and 
L — 2 sites back a further generation, in terms of the 
previous iteration |<!> i-4 ): 

\m*?n) = Y,< n \ m0 *p 0n > 
p 

1 

+ OS?" 1 In) 

r 

+ ^C>11 > f" 2 ln) (4) 

s 

where m,n = 0,1 and the L — 4 superscripts have been 
dropped for the sake of clarity. It is now possible to cast 
the Hamiltonian in a corresponding form. This involves 
some tedious algebra. However, bearing in mind that 
the states for different N are orthogonal as arc the sets 
of eigenstates $p ~ 2 and <E>g ~ 4 of the Hamiltonian at the 2 
previous iterations, the final Hamiltonian may be written 
as the following 4x4 block form: 
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For each particle number with a set of basis states, this 
block matrix can be used to generate the elements of the 
full matrix. A ground state can be calculated for each 
of these particle numbers. The ground state lowest in 
energy is the system ground state (this is how \x controls 
the ground state particle density). 



3. Reducing the Number of Basis States 

The purpose of reformulating the basis states and in 
turn the Hamiltonian in this manner is to enable an ap- 
proximation to be introduced which keeps the dimen- 
sion of the Hilbert space roughly constant as sites are 
added. During each iteration a proportion of the basis 
states must be thrown away according to some system- 
atic method. This is necessary to keep the calculation 
to a computationally manageable size. Within the tight- 
binding framework it is the only approximation in our 
method. 

There are several possible schemes which could be 
used. A criterion is required that produces the smallest 
error on the next iteration ground state as it is the prop- 



erties of the ground state which are of interest. Naively, 
the lowest energy states could be kept. More sophisti- 
cated approaches would determine which states make the 
largest contribution to the next iteration ground state. 
Whichever method is adopted, some justification will be 
required as will the determination of the limits of its ac- 
curacy. 



The simplest method to implement is to throw away 
the states of highest energy, so this will be adopted 
initially. The diagonalization routines automatically 
sort eigenstates according to their eigenenergies making 
the procedure relatively straightforward. Thus during 
each step, after diagonalization but before extending the 
chain, the highest energy states are discarded. This is 
achieved by setting an energy cutoff halfway between the 
Mth and (M + l)th eigenvalue with the same occupa- 
tion number as the system ground state. This is demon- 
strated in fig. [21 Then all states with energy higher than 
the cutoff are removed. The value of M can be changed 
to control the accuracy where higher accuracy of course 
entails larger processing time and memory requirements. 
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FIG. 2: This diagram illustrates the simple energy cutoff pro- 
cedure for reducing the Hilbert space. It shows an example of 
basis states scattered according to energy and particle num- 
ber. Firstly, the ground state is located. Then the excited 
states with the same occupation number are counted in order 
to set the cutoff energy as the midpoint between the Mth and 
(M + l)th state (in this case M = 9). All states below the 
cutoff are kept and the rest are discarded. 



C. Measuring the Localization Length 

The aim of this new method is to understand the ef- 
fect of varying system parameters, in particular the in- 
teraction strength [/, on electron localization. For non- 
interacting systems and assuming the wavefunction has 
an exponentially decaying envelope, the degree of local- 
ization can be characterized by the localization length A. 
This quantity is related to the transport and conductivity 
of the system. 

In the non-interacting case a number of methods 
have been implemented for extracting the localization 
length^. However, most methods of gauging localiza- 
tion cannot be simply carried across into the many-body 
case. The sensitivity to boundary conditions (BCs)i2ii£ 
does not suffer from this problem. 

This method works by bringing the ends of the chain 
together to form a ring, which is equivalent to changing 
the boundary conditions from open to periodic. In fact, 
in so doing it is possible to introduce a complex phase 
factor. The essence of the method is to calculate the 
change in ground state energy as the boundary condi- 
tions are twisted, i.e. as the phase factor is varied. For 
a localized state very little change should be observed as 
the wave amplitude has decayed to zero. However, an ex- 
tended state should experience a substantial change since 
the wave amplitude has not decayed off. 

The simplest way to implement this approach is to use 
the phase sensitivity D. This is defined as the differ- 
ence in ground state energy Eq between the system with 
periodic BCs and the system with anti-periodic BCs: 

D = Eq (periodic) — Eq (anti-periodic). (5) 



The length dependence of the phase sensitivity can be 
used to define A: 

Doce^. (6) 

As the new recursive method works by extending a chain 
with open boundary conditions, it was necessary to im- 
plement the phase sensitivity perturbatively. This was 
done both as an analytical perturbation and as a numer- 
ical perturbation. The former actually reduces to calcu- 
lating certain elements of the density matrix. 

1. Analytical Perturbation 

To implement the phase sensitivity as an analytical 
perturbation, consider the first order energy shift of the 
ground state |4> ) 

5E = <$ |A#|*o), (?) 

where SH = ±V(c\c L + c^CjJ + U(c\c 1 )(c' L c L ) contains 
the hopping and interaction terms now connecting the 
two ends of the chain. Calculating the effect of SH mo- 
tivates adding sites simultaneously to both ends of the 
chain. Given that |$o) is a linear combination of basis 
states (sec cqn. the effect of SH on |$o) is 

5H\^- N ) = 

± y(-i)^- i ^{6 0j |o<i> J i " 2 ^- 1 i) 

3 

+ CQ^f-^O)} 

+ Uj2dok\l^- 2 ' N - 2 l). (8) 
k 

where (— l) Si_1 is a phase factor due to electron anti- 
symmetry arising out of the occupancy of the inner L — 2 
sites. Substituting this into the expression for SEq gives 

SE Q = ±2V(-1) SL - 1 J2 b °i c °i + U J2^ 2 - ( 9 ) 

3 k 

The phase sensitivity is the difference between periodic 
and anti-periodic energy shifts, thus the last interaction 
term cancels yielding 

D = w(-iy^J2 b °3 c °3- ( 10 ) 

3 

The information required to calculate the factor (— \) SL - 1 
is unavailable. Fortunately the degree of localization can 
be calculated by using D 2 instead. Alternatively, this 
problem can also be avoided by choosing a different order 
for applying creation operators. 

The "scalar product" quantity in the expression for D 
corresponds to calculating the off-diagonal element of the 
reduced density matrix, 

P<i-o|,|0-i> = P(o-i|,|i-o) = yifrojCpj- (11) 

3 
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This is intuitively unsurprising because the quantity of 
interest is the probability that given an electron is placed 
at one end of the chain, an electron comes out the other 
end. Furthermore, only information about the ends of 
the chain is available (and required), so a reduced density 
matrix is used that sums over the redundant middle part 
of the chain. 



2. Numerical Perturbation 

The second way the phase sensitivity to boundary con- 
ditions was implemented uses a numerical perturbation. 
That is, the calculation proceeds as normal with open 
boundary conditions. At each step, working with the nor- 
mal basis states two additional Hamiltonians arc formed 
corresponding to periodic and anti-periodic BCs. These 
are then solved and the ground state for each type of BC 
is found. The phase sensitivity is then easily calculated. 
The basis states generated by the open BC Hamiltonian 
are kept as normal, but the states from the other two 
BCs are discarded. 

Using this approach means performing extra diagonal- 
izations, so it is computationally time consuming. Hence 
it was used to numerically verify the legitimacy of the 
analytical perturbation (i.e. the off-diagonal element of 
the reduced density matrix). 

The Hamiltonian for (anti-)periodic BCs is identical to 
the open BCs Hamiltonian B 2|) but with a few extra 
terms. These additional terms are those calculated in 
the analytical perturbation for each basis state JSJ and 
appear in the diagonal matrix elements of the Hamilto- 
nians. 



D. Model Properties 

The Hamiltonian models used with this new method 
possess some shared properties. These will be defined 
in this section, again in terms of the single chain model, 
and suggest some consistency tests which can be used to 
provide limited justification for the method. 



1. Definition of Energy Gaps 

Several definitions of the energy gap exist in many- 
body systems. The numerical method under develop- 
ment works within the grand canonical (GC) framework 
so eigenvalues are grand canonical energies. Thus to con- 
vert to canonical (C) energies the following relation must 
be used 



Various energy gaps are defined in i|13|) 

AE ph = E^-EoiN) 

AE+ = E (N + 1)-E (N) 

AE_ = E {N) - E (N - 1) 

AE = E a (N + 1) + E Q (N - 1) 



(13a) 
(13b) 
(13c) 
2£*,(iV).(13d) 



E (N) ~ E (N) — Nfi. 



(12) 



where AE'ph (I13afl is the difference between the ground 
state and the 1st excited state, AE + and AE_ (|13bl & 
I13c|) are the energies to add or remove and electron from 
the system, and AE is the difference of AE + and AE_ 
(|13d(l . When transferring to canonical energies, AE + and 
A_E_ will shift by a constant =F^i. The non-interacting 
limit yields some predictions which simulations should 
reproduce. 



2. Length Dependence 

In the single body case there is one state per site and 
the bandwidth is constant so the density of states is pro- 
portional to the number of sites L. This means on av- 
erage AEph oc at the same position within the band. 
The canonical and grand canonical versions are identi- 
cal as AEph is a difference between states with identical 
particle number. 

The canonical energy required to add a particle should 
be a little larger than the chemical potential. In the 
ground state, the highest occupied state will be the first 
state below the chemical potential. In order to add a 
particle the energy of the first state above the chemical 
potential is required. Thus the canonical energy is the 
chemical potential plus a contribution with a length de- 
pendence of j-. The maximum of this contribution is 
AEph corresponding to the case in which, before adding 
the electron, the highest occupied state energy was pre- 
cisely the chemical potential energy. When transferring 
into the grand canonical scheme AE + reduces to the cx 
contribution. In the middle of the band, where fi = 0, 
the canonical and grand canonical energies are identical. 



3. Consistency Test 

The energy gap definitions also give a consistency test 
that works in the presence of disorder, but unfortunately 
does not work in the presence of Coulomb interactions. 
Consider the energy change when one electron is added 
to the ground state. It can readily be seen that a relation 
can be found for the difference between ground state en- 
ergy of the system Eq (N) and the ground state with an 
additional particle Eq(N + 1) 

E Q (N + 1) - E (N) = E^N) - E (N - 1), (14) 

where E\ (N) is the first excited state with the same par- 
ticle number as the system ground state. Figure [3] shows 
the four states referred to in the relation. The right-hand 



6 



O O o 

oooo 

(5 O 

o 6 • • 
o • o • 



N-l N N+l 

^(N-l) EtfN) E,(N) EtfN+l) 

ground ground first ground 
state state excited state 
state 

FIG. 3: Illustrates the consistency test relation between 
ground state energies. The filled dots are occupied single 
particle states. 

side can be substituted into the definition of AE, then af- 
ter canceling AE reduces to AE p h- Therefore the above 
relation is equivalent to 

AE = AE ph . (15) 

The accuracy to which these equivalence relations will 
be satisfied depends on the number of basis states re- 
tained in each generation. If all states are kept this 
should be an exact formula. 



4- Particle-Hole Symmetry 

When a band is nearly full it is often more instruc- 
tive to consider the system in terms of holes rather than 
particles. Even in the case with disorder and interac- 
tions there should be a symmetry between the behavior 
of holes in the top of the band and the behavior of elec- 
trons in the bottom of the band. The symmetry relation 
between the two can be found by rewriting the Hamil- 
tonian in terms of holes and then comparing with the 
original electron Hamiltonian. 

The Hamiltonian for the system given in terms of elec- 
tron occupancy is (for open boundary conditions): 

L L-l 

i=l i=l 

L-l L 

+ u ^2(4 d i)(4+A+i) - ME e ta- ( i6 ) 

t=l i=l 



To rewrite this for holes it is necessary to define suitable 
creation and annihilation operators for holes, b\ and S, 
respectively. The equivalence relations to particle oper- 
ators are: b\ = c i (— l) 1 and "b i = c\{— 1)\ Substituting 
these into the electron Hamiltonian and rearranging gives 
a similar form of Hamiltonian to the original 

L L-l 

i=i i=i 

L-l L 

+u E( s &)$+A+i) + G" - 2U ) E s & 

i=l i=l 

L 

+ ^2ei + (L-l)U-Lfi+ Ubl^ + Ub[b L {17) 

i=l 

The last two terms arise from using open boundary con- 
ditions. In order to correct for these terms, from now 
on the original Hamiltonian will be adjusted by adding 
+^c\c x + ^c' L c L . In the algorithm these terms are read- 
ily introduced but must only be applied to the end sites, 
i.e. they must be removed when extending the chain 
length. Once added, the equivalent Hamiltonian in terms 
of holes then becomes 

L L-l 

i=l i=l 
L-l 

1=1 

+ f ^l 

L L 

- 2io E % + E e > + L{ v - ( 18 ) 

»=1 i=l 

This symmetry means that an exact correspondence of 
eigenvectors is expected to be observed at the top and 
bottom of the band such that 

H(e i =e' i ,n = -iJ,',U = U') = 

H(e l = -e\, n = n' + 2U', U = U 1 ) 

L 

+ ^ £ ' t +i([/' + M '), (19) 

i=i 

with a shift of eigenenergies according to the last two 
terms. In fact, the + Yli=i + + terms represent 
the energy of the full electron band which can be easily 
seen from the electron Hamiltonian. The significance is 
that the hole picture starts from the top of the band and 
works downward. 

Moreover, because the eigenvectors are identical, the 
reduced density matrix method for calculating localiza- 
tion length will give exactly the same results either top 
or bottom of the band. This is because pto—i\,\i—o) = 
P(i- o|,|o - i) from eqn.^] and so can be used as a consis- 
tency test. 
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This symmetry may be extended from a single sample 
to the average case by noting that the probability distri- 
bution for e, is a symmetric function about the origin. 
This has two implications. Firstly, the Yli=i £ i t erm wm 
tend to zero as the number of sites is increased. And sec- 
ondly, the sign of the — 2<=i H^i £ i term can be flipped. 
This simplifies the equivalence relation to 

H(fi = -fJ,',U = U') = 

H(p = // + 2U',U = U') + L(U' + //). (20) 

for an ensemble average. 

The particle-hole symmetry also gives the condition 
for half- filling: \j! = U' . This can be seen by setting 
— (jf = // + 2U', so that in both eans. I19landl20lthc same 
value for fx is specified. 

E. Computational Implementation 



Solve the 2 site system 



V 

"W Remove the high energy states 



f 

Extend the length of the chain 



f 

for each particle number: 

[ Calculate the Hamiltonian matrix elements 
Diagonalize the Hamiltonian 1 



\ 

Extract physical quantities of interest 



FIG. 4: A flow diagram showing the central procedure of the 
algorithm. This whole procedure is repeated many times with 
different disorder realizations. 

The simulation was written in CH — h, making use of 
the object orientated facilities. Because states need to 
be added and removed from a set of states it was natural 
to use linked lists, with each link holding data for one 
state and the entire list representing a set of basis states 
or eigenstates. This results in making the code for the 
central algorithm less cumbersome. The structure of the 
central algorithm is shown in fig. 0] 

To extract sensible data localization quantities must 
be averaged over many systems. Conventional practice 
is to perform a geometric average which is achieved by 
averaging the logarithm of the phase sensitivity. Then 
least-square fits were carried out to extract the local- 
ization length over a minimum of 10 sites. In addition 



other quantities were recorded such as the particle den- 
sity, ground state energy and energy gaps. 



1. Removing States 



energy " 











































- ~~~~~ keep this 
state 

no states 
below 
cutoff 



N-4 N-3 N-2 N-1 N \ N+1 N+2 N+3 N+4 N+5 P """ 



J 
ground 

Mole 

FIG. 5: This diagram illustrates how the simple energy cutoff 
scheme may occasionally cause difficulties. All states with 
iV + 3 particles will be removed, however some states with N+ 
4 particles will be retained. To avoid considerable overhead 
in coding, the lowest energy state for N + 3 will be kept. 

One minor modification to the above procedure was 
made. It is simplest to code this new numerical method 
assuming that states exist for every particle number be- 
tween two limits (in a range roughly centered around the 
ground state particle density). For example in fig. El af- 
ter the Hilbert space reduction, each value of N between 
N — 2 and N + 4 inclusive have states remaining. With 
the reduction scheme outlined above a scenario occasion- 
ally arises whereby this assumption would be rendered 
invalid. Figure demonstrates that it is possible for one 
value of N, in the fig. N + 3, to have all states removed 
yet neighboring particle numbers to have states left. To 
cater for this rare event the code would become unneces- 
sarily complicated. The problem can be easily overcome 
in such cases by retaining the lowest state even though 
it is above the energy cutoff. 

In addition to this, the ground states corresponding to 
A^ + 1 and A^ — 1 particle numbers were always retained 
at each iteration. This ensures the various energy gaps 
can be calculated even in the rare event when the N — I 
or N + 1 ground state lies above the cutoff energy. 



III. THE SINGLE CHAIN MODEL 

This method was first applied to the Hamiltonian l|TJl 
plus the two correction terms to ensure particle-hole sym- 
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mctry: 
H 



E 



L-l 



c\c i e l 



•i+l 



L-l 



U 



-c c 



1°1 



7T C L C L 



L 

-^E £ ^ 

i=l 



(21) 



This is the conventional ID Anderson model with nearest 
neighbor interactions. In this paper hopping will always 
be set to V — 1 (hence defining the energy scale for both 
U and W). The following two subsections outline some 
useful features already known about this model in two 
limits: no interactions (U = 0) and no disorder (W = 0). 



A. Non-Interacting Behavior 

A single chain with one orbital per site has one band 
centered on zero with bandwidth AV. Including disor- 
der will blur the edges, effectively widening the band. 
The localization properties of a one-dimensional non- 
interacting chain are well established. For any amount of 
disorder all eigenstates are localized. The dependence of 
localization length on disorder is usually quoted &sH 



A- 1 = 



W 



24(4y 2 - n 2 ) ' 



(22) 



This is only valid for small disorder. Note that the local- 
ization length diverges in the clean limit. Therefore, an 
important test for the new recursive method is to repro- 
duce this behavior. However, care is required in making 
the correct comparison: how does this dependence carry 
across from the single-particle case to the many-particle 
case? A simple test program was constructed using ex- 
act diagonalization to calculate the phase sensitivity to 
boundary conditions for both the single state in the cen- 
terer of the band and the phase sensitivity for the grand 
canonical energy. Good agreement was found. 



At half-filling there arc two limiting forms of the 
ground state with a crossover regime. For large repulsive 
interactions a charge density wave (CDW) is observed 
(i.e. alternate sites are occupied). For attractive interac- 
tions, electrons tend to cluster together. For open bound- 
ary conditions (with the correction for particle-hole sym- 
metry) the transition region occurs between U = —0.8 
and U = -2. 




Particle Number (N) 

FIG. 6: Results from a (clean) short chain of 10 sites demon- 
strating phase separation for U < — 2. For each particle num- 
ber the ground state energy is plotted. The chemical poten- 
tial is set to give half-filling as the overall ground state (i.e. 
H = U). Plotted energies are grand canonical. 

For attractive interaction below U = —2, it is impos- 
sible to maintain half-filling within the grand canonical 
scheme. The ground state is a completely empty or com- 
pletely full band (i.e. it is unstable to phase separation) 
as can be seen in figure In fact, as the U = — 2 limit 
is reached from above, the ground state energy tends to- 
ward being independent of particle number N. 

In contrast, for increasing repulsive interactions, above 
U = 2 a charge gap opens upi^i^. In other words, the 
CDW above this point corresponds to a Mott insulator 
state. For short chains it is not possible to pinpoint where 
this gap begins. 



B. Clean Phase Space 



Previous Work 



The second limit to be outlined is the zero disorder 
phase space. This is understood because without ran- 
domness the present model can be mapped to a XXZ spin 
chain model and solved exactly for half-fillingiii^iSiii. 

A program was constructed to perform exact diago- 
nalizations on short chains in order to compare results. 
This was necessary because the method under develop- 
ment doesn't use the particle occupation basis. The com- 
putational limit is about 10 sites, but nevertheless gives 
valuable insight into the nature of the ground state for 
different regions of phase space. 



Perhaps the most substantial work is that of Gia- 
marchi and Schulz 1 ™, although they acknowledge an ear- 
lier papei— . For the present one-dimensional model 
repulsive interactions increase localization because the 
CDW is pinned by the disorder. In contrast, attractive 
interactions decrease localization. In fact, a delocalizcd 
phase is predicted for sufficiently attractive interactions. 
Giamarchi and Schulz develop a fc-space Renormalization 
Group approach to study the transition. The existence 
of this transition is ascribed to competition between dis- 
order and superconducting fluctuations^. 
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Numerical work has sought to verify these predictions 
and in particular to map out the dclocalizcd regime. One 
paper— performs exact diagonalizations on small systems 
(up to 22 sites). The results are consistent with the 
expected delocalized phase, although because the chain 
length is so small they cannot exclude the possibility that 
the localization length is very large. 

The first DMRG studyii focused on the effect of dis- 
order on the Mott state. The authors conclude that even 
weak disorder destroys the charge gap and long-range or- 
der associated with the CDW state (although the nature 
of elementary excitations remain unchanged) . 

The most extensive work has been conducted by 
Schmittcckert et al. applying DMRG to both the inter- 
acting Anderson model and to the related problem of 
persistent currents in mesoscopic ringsLi^SiiSS. The first 
study examining Anderson localization^ was on chains 
extending up to 60 lattice sites. The degree of localiza- 
tion was measured by the phase sensitivity to boundary 
conditions. Two regimes were found: a localized phase, 
U > — 1, and delocalized regime, U < — 1, consistent 
with work already mentioned. In fact, it was found re- 
pulsive interactions increase localization. After consid- 
erable numerical effort a phase diagram was produced 
showing where the two regions lie in disorder-interaction 
space. The authors believe the earlier attempt— based on 
an RG procedure over estimates the delocalized regime 
by a factor of 4. Other authors, Romcr— and Schus- 
ter ct al— , have mapped out an extended regime for the 
same model but with the Aubry- Andre quasi-periodic po- 
tential. However its shape in disorder-interaction phase 
space takes on a different form. 

In two more recent papersS^ Schmittcckert et al. 
showed that important physics is washed out in the pro- 
cess of ensemble averaging used in ascertaining numerical 
data. They examined the chaotic region between a chain 
characterized as an Anderson insulator and characterized 
as a Mott insulator corresponding to the strongly disor- 
dered and the strongly correlated limits respectively. 

Schmittcckert ct al. showed that in this transition re- 
gion there are sharp charge reorganizations causing a dra- 
matic increase in phase sensitivity (delocalization) of up 
to 4 orders of magnitude. It is important to note that 
the position in parameter space where these reorganiza- 
tions occur is sample dependent such that on averaging 
over many samples only a slight delocalization effect can 
be observed. This is found for large disorder W = 7, 9 
and for small repulsive interactions U < V. Typically 
two or three distinct charge reorganizations may occur 
per sample as it is moved between a Mott insulator and 
an Anderson insulator state. Attractive interactions fa- 
vor a more inhomogeneous charge density (forming clus- 
ters) and repulsive interactions favor more homogeneous 
charge density (CDW) as expected. Of course, this ef- 
fect will be most pronounced at half filling since for lower 
densities electrons can avoid each other spatially. 

This phenomenon was explained further in terms of 
avoided level crossings^ of the ground state and first- 



excited state. However, to be critical these results are 
based on small system sizes (20 sites) and may well be 
due to finite-size effects. 



D. Determining Limits and Accuracy 

As appropriate for any new method, the extent of its 
validity should be tested before producing results. There 
are three points which should be established: 

• The program is working properly. 

• The two methods for gauging the degree of local- 
ization yield consistent results. 

• The approximation for eliminating basis states is 
controlled. 

The last of these will be the most involved. 



1. Determining Numerical Limits 
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FIG. 7: The graph is a plot of the average of ^ ln(D 2 ) as a 
function of chain length L where D is the phase sensitivity 
(11 01) . Averages are geometric, that is the mean of ln(_D 2 ) is 
found. The chain length is allowed to extend until numeri- 
cal precision is lost. The /i = case was stopped at 20,000 
sites as it had saturated long before. System parameters are 
W = 2, U — and the energy cutoff is set using M = 10 
(corresponding to a total of approximately 160 basis states 
per iteration). The averages were taken over 2000 systems. 



The first simulations aimed to determine the broad 
numerical limits of this new method. Figure shows 
some typical results where the chain length has been al- 
lowed to extend as far as possible. A range of values of \i 
were used, spread across the band (between fj, = — 2 and 
fx = 2 in the non-interacting case). The lines pair up, 
with (i = —fx 1 and ft = +fx' giving near identical results 
and thus demonstrating the anticipated symmetry in the 
band. 
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Straight lines indicate exponential localization, which 
is seen near the band edge. However, in the center 
of the band different behavior is observed: the curves 
show some form of decay which saturates at large chain 
lengths. In fact, for the fi = case, the curve had clearly 
saturated and so was stopped at 20 000 sites. 

Apart from the \x = case, all simulations contin- 
ued until numerical precision limits were reached. This 
limit is encountered when calculating the phase sensitiv- 
ity. Off-diagonal elements of the reduced density matrix 
correspond to "scalar product" type quantities ljTT|) . Us- 
ing this as an analogy, when the two vectors become al- 
most perpendicular the product becomes very small. In 
this limit, numerical rounding dominates over the physics 
rendering any results based on this regime meaningless. 
Arising out of this, a criterion was devised to automati- 
cally halt simulations before this numerically inaccurate 
region is reached. Taking up the analogy again, when 
the scalar product divided by the norm of the two vec- 
tors is comparable to the floating-point precision then 
only "noise" is being calculated. This condition gives 
the upper length limit for linear fits which determine the 
inverse localization length. A lower limit was also set 
in which typically the first 20% of sites were ignored to 
allow the simulation to "settle down" . 



2. Particle-Hole Symmetry Test 

The particle-hole symmetry test can be verified by 
looking at two chains using the same random distribu- 
tion of site energies, but with one the negative of the 
other (see ean. lliJf) . On doing so identical results are ob- 
tained. Thus in the non-interacting case the electron-hole 
consistency test is convincingly satisfied. 

Applying a small electron interaction force U = 0.1 
causes the paired lines to split. This is also expected. 



Once the + ^■c\c 1 + ^c\c L correction terms are included, 
the two lines then collapse on top of each other again. 



3. Comparison of the two Phase Sensitivity Methods 

As explained earlier in section lTl CI the reduced density 
matrix method of determining the extent of localization 
could be verified by calculating the phase sensitivity as a 
numerical perturbation. It was found that the two meth- 
ods are in good agreement (fig.|Hl)- Given the significant 
computational processing time required for the numerical 
perturbation method, normally only the reduced density 
matrix method will be used. 



4- Reducing the Number of Basis States (Revisited) 

The initial results (fig. 0) show that the method fails 
in the middle of the band - exponential decay is not ob- 
served in the non-interacting case. This must be due to 
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FIG. 8: This graph shows the inverse localization length (ob- 
tained from a linear fit) against a range of values of the chem- 
ical potential fi. The two lines correspond to the two methods 
of determining localization and they show very good agree- 
ment. System parameters are W — 2, U = with chains 
allowed to extend to a maximum of 500 sites. The averages 
were taken over 2500 systems and the energy cutoff set by 
M = 10 (corresponding to a total of approximately 160 basis 
states per iteration). The error bars represent one standard 
deviation. 
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FIG. 9: The graph is a plot of the average of | ln(D 2 ) against 
chain length L. The parameters are the same as fig. Q but 
restricted to the middle of the band (fj, — 0). The only differ- 
ence between the three curves is the procedure for discarding 
basis states. 



the Hilbert Space reduction criterion, as it is the only ap- 
proximation in the method. A simple variant on the orig- 
inal procedure was tried: the total number of states for 
all particle numbers was fixed rather than using a fixed 
number for the ground state particle number only. This 
simple modification induced a significant change in the 
results. Although the decay was still non-exponential, 
this change clearly yields an improvement toward the ex- 
pected behavior (fig. EJ. 

Presumably the key difference between the methods is 
that using states from all particle numbers results in an 
energy cutoff which fluctuates less. The natural criterion 
to try next is a fixed energy cutoff. This can be done by 
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averaging the value of the cutoff using the fixed number 
of states method. Then a fit of the cutoff as a function 
of chain length could be used as a fixed energy cutoff. 
It turns out that it is not possible to do this as an ab- 
solute cutoff because the ground state energy fluctuates 
too much. However it can successfully be done as a fixed 
energy cutoff relative to the ground state. When imple- 
mented exponential decay is observed in the middle of 
the band (fig.®. 

We conjecture that this dramatic improvement, result- 
ing from an apparently innocuous change in cutoff meth- 
ods, can be explained in terms of energy level statistics. 
Consider the middle of the band with no interactions: 
electrons should be localized, with states obeying Pois- 
son statistics. One may envisage the system accidentally 
encountering a higher density of low lying energy states. 
According to the original method, the energy cutoff is 
correspondingly lower. Thinking in terms of energy level 
repulsion, this would result in a release of "pressure" as 
a larger number of states are removed. The opposite 
scenario in which states accidentally spread wider than 
average may also be considered. In this case, the cutoff 
has a smaller effect than normal. The combined effect is 
to reduce fluctuations, causing the system to bear more 
resemblance to a Wigner distribution. In others words 
the system tends toward derealization, consistent with 
the data on fig. [5] 

The second cutoff method implemented worked by fix- 
ing the cutoff using states across all particle numbers. 
According to the picture just outlined, the same effect of 
dampening fluctuations should still be present, although 
less severe because using a greater number of states re- 
duces fluctuations of the cutoff energy. This can also be 
observed in figure El where the delocalizing effect is not 
so strong. 

The third procedure for discarding states used a fixed 
energy cutoff, which is completely uncorrelated to the 
density of low lying states. Therefore the delocalizing 
effect is completely absent. 

E. Comparison with Non-Interacting Results 

The new fixed energy cutoff may be used to pro- 
vide further verification by checking that when interac- 
tions are turned off non-interacting results can be repro- 
duced. This is particularly important as some well es- 
tablished methods applied to the two-interacting particle 
model can fail in this respect (e.g. the Transfer Matrix 
Methoc&2&). 

1. Dependence on Disorder 

As a test of the accuracy of the method, figure ITUI was 
produced. The calculated results should correspond to 
the known result (|22() . The new method seems to give a 
dependence on W greater than W 2 . Values are at least 
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FIG. 10: Graph showing the dependence of inverse localiza- 
tion length upon the disorder when interactions are turned 
off. Results from the new method should correspond to known 
result for the middle of the band. Systems were allowed to 
extend up to 1000 sites, retaining an average of 480 basis 
states per iteration. Averages were taken over 1000 disorder 
realizations. 
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FIG. 11: Graph showing single sample results for low disorder 
(W = 0.2). The new method fails in this limit. The only 
difference between the lines is a slight variation in the energy 
cutoff. 



of the correct order of magnitude. It should be noted, 
however, that in the clean limit and for low disorder the 
new method fails to produce meaningful results. Fig- 
ure ^2 shows a set of results corresponding to the same 
system but with the energy cutoff slightly varied. The 
change in behavior is quite dramatic and dominated by 
large oscillations. These effects do not occur for stronger 
disorder. 
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FIG. 12: Graph showing the dependence of inverse localiza- 
tion length on the average number of basis states retained per 
iteration for the non-interacting case with disorder W — 2 
(left) and W — 5 (right). The insets show this quantity plot- 
ted against the actual energy cutoff used. Such plots can be 
used to test for convergence. Data was averaged over 100 
systems with chains extending up to 1000 sites. 



2. Convergence 

Figure E| is used to determine whether the method 
converges to the known value of the localization length 
in the middle of the band. For W = 2, cqn. [221 im- 
plies A -1 ps 0.038. Hence for the calculation using the 
largest matrices the localization length is over estimated 
by about a factor of 4. As anticipated from fig. ^3 con- 
vergence is much better for W — 5. In this case, ean. 1221 
gives A -1 ps 0.24. Figure IT2"1 shows the new method con- 
verging at a value close to this result. One could con- 
sider proceeding by just examining interaction effects for 
W > 4. However, interesting physics is expected when 
disorder and interactions are of similar strengt. Note 
that the standard deviation, W/ \/\2, is a more satisfac- 
tory measure of disorder, so that this condition is fulfilled 
when W = a/12 and U = ±1. 




FIG. 13: Graphs showing the dependence of the charge gap 
AE on chain length. The black lines are results from the new 
method under development and for comparison the 2nd line 
represents exact results from a single-body calculation. The 
left-hand side shows the clean case and the right-hand side is 
the average over 1000 systems with disorder W = 2. The new 
fixed energy cutoff was used to remove basis states (with an 
average of 700 per iteration). 



IV. RESULTS 

Despite the unanswered questions, results were suc- 
cessfully obtained using the fixed energy cutoff procedure 
for eliminating states. The value for the cutoff was de- 
termined by first using the fixed number procedure for 
a small number of systems. When the average number 
of basis states is given, it refers to the number of basis 
states used with the fixed number procedure in order to 
determine the fixed energy cutoff. In addition note that 
previous work by other authors focuses on the case of 
half-filling, so results presented in this section also exam- 
ine this particle density. 



3. Energy Gaps 

Earlier, in section III D II it was noted that in the non- 
interacting limit there is a relation between two of the 
energy gap definitions (15} . This condition was found 
to be obeyed provided a large proportion of basis states 
were retained (i.e. only possible for short chain lengths). 

Also in this limit, the energy gaps should have a spe- 
cific dependence on length. This is rigorously known in 
the clean limit as oc \ for AE, AE p h and AE + . In 
fact, when averaged, systems with disorder exhibit the 
same behavior. Figure 1131 shows how results from the 
new method compares with the non-interacting results. 
In both cases the decay was only satisfactorily observed 
for short chain lengths. The charge gap decays to a value 
above zero. The oscillations arise due to the removal of 
states which can be established by varying the number of 
states retained. This saturation and oscillatory behavior 
was also observed for the other two gap definitions. 



A. Dependence on Interaction Strength 

Of central interest is the effect of electron-electron in- 
teractions on localization. FigurcHHshows the calculated 
dependency for disorder W = 2. Three different energy 
cutoffs were used corresponding to different numbers of 
retained states. 

The overall behavior is unambiguous: repulsive inter- 
actions enhance the effect of disorder whereas attractive 
interactions reduce it. For U > the inverse localization 
length has an approximately linear relationship to inter- 
action strength. Below U = — 1 the localization length 
diverges. This apparently extended regime is anticipated 
from previous work (section ITTT"0|) . 

Note that due to the "flattening" effect (fig. EJ) toward 
the phase separation the many-body density of states 
rises rapidly with energy. For a fixed energy cutoff this 
means that more states are retained over a greater range 
of particle numbers. This reduces computational perfor- 
mance and so the region —2 < U < —1.8 has not been 
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FIG. 14: The dependence of localization length on interaction 
strength. The three lines correspond to different energy cutoff 
values. Each line is averaged over 1000 systems which are 
allowed to extend to a maximum of 1000 lattice sites. Disorder 
W = 2. 
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explored. Data was obtained down to about U = —1.8 
and is displayed in fig. IT"!! 




1000 2000 3000 4000 

Average Number of Basis Stales 



FIG. 15: Graph showing the dependence of inverse localiza- 
tion length on the average number of retained basis states 
in order to test convergence in the delocalized regime. The 
value of interaction strength chosen was U = — 1.4. The inset 
shows the dependence on the fixed cutoff energy. The data is 
averaged over 100 systems with chains allowed to extend to a 
maximum of 1000 sites. Disorder W = 2. 

To explore the convergence, one point in the middle 
of the delocalized regime was chosen and its convergence 
properties were explored. Figure ITol demonstrates that, 
within the computational limits, the extended regime is 
robust. 



B. Disorder-Interaction Phase Space 

Having confirmed the existence of a delocalized regime 
for attractive interactions, it was natural to attempt to 
map out the extent of this region. This was done in 
disorder-interaction phase space as plotted in figEH The 
area marked as delocalized corresponds to systems with a 



FIG. 16: Disorder-interaction phase space plot for the sin- 
gle chain model at half-filling. The spectrum represents the 
degree of localization. The lowest interval corresponds to a 
localization length greater than 1000 sites. This contour plot 
was produced using over 1300 points. Each point was aver- 
aged over 250 systems in which chains were allowed to extend 
to 2000 sites and approximately 240 basis states were retained 
per iteration. Data for W < 0.6 is not shown because the 
method is unreliable for low disorder as discussed earlier. 

localization length greater than 1000 sites. If the number 
of systems averaged over were increased then this crite- 
rion could be made more stringent. With this definition 
it can be seen that the dclocalizing effect does not extend 
beyond W = 2.5. The extended region appears to cross 
the non-interacting line for lower disorder. However, the 
method is unable to produce meaningful results for low 
disorder. The phase sensitivity tends to be dominated by 
oscillations which appear as the clean limit is approached 
(see fie lllf) . The (unreliable) results which were obtained 
indicate that the proximity of the derealization in fact 
decreases for lower disorder. Hence this remains an open 
question. 

C. Determination of Exponents 

It is tempting to describe the change from a localized 
system to a delocalized system in the language of second- 
order phase transitions. Although, there is no a priori 
reason for believing the transition can be described in 
this way, it is clear that a description as a first-order 
transition is ruled out by the absence of any discontinu- 
ities. Working on this assumption, the obvious quantity 
to calculate is the exponent defined as 

A" 1 = A(U - U c ) v U > U c , (23) 

where A is a coefficient, v is the exponent and U c is the 
critical interaction strength where the transition occurs. 



14 




-1.2 -1 -0.8 -0.6 
Interaction [U] 



FIG. 17: This graph is the same as the black line in fig. H4l but 
focusing on the transition region. The fitted line corresponds 
to U c — —1.375 and v = 3.0 and shows good agreement. 



Two factors make determining this exponent difficult: 
Firstly, data can only be used on one side of the tran- 
sition. Secondly, the precise value of U c is unknown and 
no scaling analysis is available. Consequently, fitting the 
data shown on fig.^Jusing a log-log plot proved unsatis- 
factory. Likewise, attempting to fit all 3 parameters (A, 
U c and v) simultaneously was also uncontrolled. There- 
fore, an estimate for the critical interaction strength was 
visually estimated and then the other two parameters 
could be fitted. Although, this procedure is less than 
ideal, it appeared to be the most satisfactory approach. 
The estimated range for U c is —1.45 < U c < —1.3. Us- 
ing a number of values within this range the two other 
parameters (A and v) were fitted, giving v ss 3.0 ± 0.4. 
This exponent has been plotted with the central value of 
U c on fig. El The fit is not perfect, the deviation may 
either reflect the absence of a known scaling analysis or 
that the assumption of a second-order phase transition is 
incorrect. 

However, it is at least clear that away from the transi- 
tion there is an approximately linear relationship between 



inverse localization length and the interaction strength. 
Further, by inspecting fig. 1141 the exponent must there- 
fore be greater than 1 in the vicinity of the transition. 
The determined value for v is consistent with that obser- 
vation. 

It would be desirable to perform a similar procedure 
for other values of disorder and to determine the disorder 
exponent at a fixed interaction strength (i.e. approach- 
ing the transition from above on fig. I16|) . However, the 
quality of data presently obtained is inadequate. Schmit- 
teckert et al^, argue that the exponent is non-universal. 



D. Summary 

We have presented a method of studying disordered 
and interacting quasi- 1-dimcnsional systems which com- 
bines aspects of the transfer matrix and DMRG ap- 
proaches. While the method works well and is able 
to study significantly larger systems than have been 
achieved hitherto, there is still room for improvement. 
In particular the strategy for reducing the Hilbert space 
and compensating for the side effects of the reduction is 
still too simplistic. It would be useful to understand why 
the method fails so dramatically for low disorder. Nev- 
ertheless the method is generalizablc to more complex 
problems such as the Hubbard model or strips of finite 
width. It could eventually be possible to combine such 
an approach with finite size scaling in order to study the 
metal-insulator transition. 

As a first application of our method we have presented 
results on spinlcss fcrmions in ID. There is qualitative 
agreement with previous work: repulsive interactions in- 
crease the effect of disorder and attractive interactions 
have the opposite effect. We have mapped out the delo- 
calized regime and found some disagreement with previ- 
ous work. According to our results, DMRG studies under 
estimate this region by a factor of 2 and an earlier study 
over estimates it by a factor of 2. We have also made a 
first estimate of the critical exponent of this transition, 
but our data is not yet sufficient to test its universality. 
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